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I. INTRODUCTION 



The Compton scattering provides a unique tool for studying hadronic structure. In distinction to hadronic form 
factors, the Compton ampUtudc probes the hadrons through a couphng of two electromagnetic currents, and thus 
provides a wealth of additional information. When the initial photon is highly virtual while the final one is real, 
one arrives at the kinematics of deeply virtual Compton scattering (DVCS), the process [1,2] which was a subject of 
intensive theoretical analysis during the last decade. Within the perturbative Quantum ChromoDynamics (QCD), it 
was established [3,1,2,4,5] that the most important contribution to the amplitude of this exclusive process is given by 
the convolution of a hard quark propagator and a nonperturbative function describing long-distance dynamics which 
is now known as a "generalized parton distribution (GPD)" [3,1,2] (see [6-8] for recent detailed reviews). The GPDs 
serve as a generalization of the ordinary (forward) parton distributions and provide much more direct and sensitive 
information on the light-front (LF) wavefunctions of the target hadron than the hadron form factors. In particular, 
the momentum of the "probed quark" in GPDs is kept fixed at longitudinal momentum fraction x, while for the form 
factor it is integrated out (due to the nonlocal current operator ■tlj{0)'y^tp{z) for the GPDs in contrast to the local 
vertex ■i/5(0)7'''(/'(0) for the form factor). Within the GPD approach, form factors are treated as just the lowest moments 
of the GPDs. Because of being generalized amplitudes, the GPDs always involve the nonvalence contributions due 
to the presence of the asymmetry between the longitudinal initial P and final P' hadron momenta characterized by 
the "skewness" parameter = (P — P')^ /P^. The kinematic region where the longitudinal momentum fraction 
X of the probed quark is greater than the skewness parameter C, {i.e. 1 > x > C,) is called "DGLAP region" (the 
evolution pattern for GPDs there is similar to the DGLAP evolution [9-11] of ordinary parton distributions) while the 
remaining part < a; < ^ of the longitudinal momenta < a: < C is called "ERBL region" (GPDs in this case have 
ERBL-type evolution [12,13] characteristic of meson distribution amplitudes). In the framework of LF dynamics, the 
DGLAP and ERBL regions have also been denoted as the valence and nonvalence regions, respectively, because the 
parton-number-changing nonvalence Fock-state contributions cannot be avoided for < .t < C while only the parton- 
number- conserving valence Fock-state contributions are needed for 1 > a; > C. Thus, it has been a great challenge to 
calculate the nonvalence contributions to the GPDs in the framework of LF quantization. 

Although many recent theoretical endeavors [14-20] have been made in describing the GPDs in terms of LF wave- 
functions, the results have not yet been satisfactory enough for practical calculations. In Refs. [15] and [16], the 
nonvalence contributions to the GPDs have been rewritten in terms of LF wavefunctions with different parton configu- 
rations. However, the representation given in Refs. [15] and [16] requires to find all the higher Fock-state wavefunctions 
while there has been relatively little progress in computing the basic wavefunctions of hadrons from first principles. 
In Refs. [17] and [18], the GPDs were expressed in terms of LF wavefunctions, but only within toy models such as 
the 't Hooft model of (1 + l)-dimensional QCD [17] and the scalar Wick-Cutkosky model [18]. While these toy model 
analyses are helpful to gain physical insight on the properties of the GPDs (especially, the time reversal invariance, 
the continuity at the crossover between the DGLAP and ERBL regions, and the sum rule constrained by the electro- 
magnetic form factor), the real (3 -|- l)-dimensional QCD motivates us to come up with a more realistic model for the 
application to the analysis of GPDs. 

In an effort toward this direction, one of us (C. Ji) presented an effective treatment of handling the nonvalence 
contributions to the GPDs of the pion [19,20] using the LF constituent quark model (LFQM), which has been phe- 
nomenologically quite successful in describing the space-like form factors for the electromagnetic and radiative decays 
of pseudoscalar and vector mesons [21-23] and the time- like weak form factors for exclusive semileptonic and rare 
decays of pseudoscalar mesons [24 26]. Our effective treatment of handling the nonvalence contributions is based on 
the covariant Bethe-Salpeter (BS) approach formulated in the LF quantization [24] which we call LFBS approach. 
It has been previously applied to the exclusive semileptonic and rare decays of pseudoscalar mesons [25,26] provid- 
ing reasonable results compared to the data. The projection of the four-dimensional two-body BS equations at the 
light-front hypersurface has also been discussed in Ref. [27] and the problem of constructing gauge invariant current 
in terms of LF bound-state wavefunctions has been handled in Ref. [28] 

In Ref. [19], an artifact of discontinuity dX x = Q between valence and nonvalence parts of the GPD had occurred 
but it was later cured in the framework of the same LFBS approach by taking into account both the vertices of the 
meson and of the gauge boson [20] . Due to the consistent treatment of the vertices, the continuity of the GPD at the 



2 



crossover point x = C, was secured. However, the value of the GPD at the crossover point vanished due to the end-point 
behavior of the LF wavefunctions of meson and gauge boson. In Ref. [20], it was noted that it may be possible to 
avoid this effect of GPDs' vanishing at the crossover point by considering contributions of the higher Fock states in 
the DVCS amplitude. In this work, we examine if the GPD value at the crossover point is indeed nonzero including 
the higher Fock state corrections. This investigation is particularly important because the recent measurement of the 
single spin asymmetry (SSA) at HERMES [29] and CLAS [30] indicate that the value of GPD at the crossover point 
does not vanish for the proton. Although the internal structure of the pion is rather different from that of the proton, 
the issue of GPDs' vanishing or not vanishing at the crossover point is common to all bound states. Thus, establishing 
the link between the nonzero value of GPDs at the crossover point and the higher Fock-state contributions is a result 
applicable for both hadrons. 

Working in the same LFBS framework, we will consider in addition to the pion LF wavefunction the contribution 
from the higher Fock components of qqg. We first follow the analysis of the higher Fock state contributions presented 
in Rcf. [18] and find that the analysis of Ref. [18] is correct only if one uses the exact solution of the BS equation. Later, 
the authors of Ref. [18] investigated the current matrix elements in the LFBS formalism discussing the replacement 
of non- wavefunction vertices [31] and applying it to GPDs [32]. Using a simple model, they obtained expressions for 
GPDs that are continuous and noted that the nonvanishing of GPDs at the crossover point is tied to higher Fock 
components. Our work develops the QCD application to the continuous and nonvanishing GPDs at the crossover 
point. We present explicit numerical results in QCD. For a model LF wavefunction we show that one shall retain all 
contributions including the ones that ideally could be absorbed into the LF wavefunction through the Bethe-Salpeter 
equation. This suggests a necessary improvement of the model wavefunctions by iterating them with the BS kernel. 

This work is organized as follows. In Section II, we go over briefly the essentials of the LF kinematics of the 
DVCS. In Section III, we review previous advances in LFBS effective treatment of the pion GPD. Then, in Section 
IV, we extend LFBS approach to include qqg Fock states and present the details of our calculations and numerical 
results including the application of the prescription given in Ref. [18]. Concluding remarks follow in Section IV. The 
pole assignment for the Cauchy integration and the comments on the organization of the numerical calculations are 
presented in Appendices A and B, respectively. 

II. LIGHT-FRONT KINEMATICS OF THE DEEPLY VIRTUAL COMPTON SCATTERING 

We begin with the kinematics of the virtual Compton scattering off the pion (see Fig.l) 

7*(5)+7r(P)^7(9')+'r(n- (1) 




FIG. 1. Handbag diagrams that give the dominant contribution to Compton scattering in the deeply virtual regime. The 
lower "soft" part consists of a hadronic matrix element which is parametrized by generalized parton distribution functions. 

The initial (final) hadron state is characterized by the momentum P {P'), and the incoming virtual (space- like) and 
outgoing real photons arc characterized by the momenta q and q', respectively. In this work, we use the LF-metric 
V ■ V = V~^V~ — V^. Defining the four momentum transfer A = P — P' , one has 

P = (^P+, M^/P+, Oj}j , P' = (^(1 - C)P+, (M^ + Ai)/(1 - C)P+, - A^) , 
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and 

A = P-P' = (CP+, (A2 + Ai)/CP+, A^) , 

where M is the mass of the pion and ( = A+/P+ is the skewness parameter describing the longitudinal momentum 
asymmetry of GPDs. The squared momentum transfer is given by 

t = A^=2P.A = -^^f_+^^ (2) 

Since > at given C, the value of —t is constrained from below: —t > — tmin = C^M'^/{1 — Q. 
As shown in Fig.l, the parton emitted by the pion has the momentum k, and the absorbed parton has the momentum 
k' — k — A. Just like in the case of space-like form factors, we may choose a frame where the incident space-like 
photon has zero plus component, g'"'" = 0, so that 

q = (0, (q^ + A^) VC-P+ + iCM^ + Ai)/(1 - 0P+,<1±) , q' = {CP^, (q± + ^±)VCP+, q± + A^) . (3) 

In DVCS, where Q'^ = -q^ » M^, and Q"^ » -t is large, C plays the role of the Bjorken variable, i.e. Q'^/{2P-q) = 
For a fixed value of —t, the allowed range of ( is given by 



0<<<y^f,l + i^-ll. (4) 



2M2 \\j {-t) 

In the leading twist, ignoring interactions at the quark-gauge-boson (photon in this case) vertex, the dominant con- 
tribution to the Compton scattering amplitude in the deeply virtual region is given by 

M^i^ = M^" + M!^" , (5) 

i.e. the sum of the s-channel amplitude Mj^" and the M-channel amplitude M^^" shown in Fig.l. The s-channel 
amplitude is given by 

Mr = -iiVee2 j ^Tr[75(/^ + m)7''(|^+ i + A + m)75(- f + ^ + m)] 



[fc2 — m? + ie] [{k + g)^ — iri^ + ie\ \{k — A)^ — m? + ie] [{P — A:)^ — m? + ie] ' 

where is the color factor and Hcov{k,P) [H'^^{k — A,P — A)] is the covariant initial [final] state meson-quark 
vertex function that satisfies the BS equation. As usual in the LFBS formalism, we assume that the covariant vertex 
function i?cov(fc) does not alter the k~ pole structure in Eq. (6). The u-channel amplitude can be easily obtained by 

M^" = M^^iq -q'). 

The effective gauge-boson vertex is depicted in Fig. 2 as a reduction from the Compton scattering amplitude of the 
virtual photon given by 

/i{4+fi + m)/ fi + m 



^ '^"^ i {q + fc)2 -m^+ie^ {k- q'f -m^+ie' ' 




FIG. 2. Effective vertex 
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In the deep inelastic limit, we may neglect ^ along with any 4-products not involving q in the trace of the numerator 
(note that for circularly polarized photons also q ■ e and q ■ e' can be neglected), i.e. 

q^ ,q ■ a a ■ b, m^, e- q,e' ■ q. (8) 

Then the amplitude (7) can be rewritten as 

_,e^ ( ^All + IZI ^ (9) 

" \{q + kf -m^ + ie {k- q'f -m? + ie) ' ^ ' 

Furthermore, / / can bo written as e • e' + ?'(T^'^e,j,eIy, whore cr'"^ = jjh''. 7'^]- The axial term can be neglected because 
it vanishes after integration over kj_. Thus, for the reduced vertex we obtain 



-ieli-e ■ e') 4 
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{q + fc)^ — w? +i£ {k — q')"^ — m? + ie 



The denominators can be further simplified by 

( I \ ^ . ^ ( ' + . (11) 

\{q + k)'^ — rm? + ie [k — q'Y — w? + is J P~^q \x — + is x — ie J 

Adding s- and u-channel amplitudes, we obtain the Compton scattering amplitude in DVCS limit as follows 

M" = = -|jeS / |t (^-i^ + P)HUt - A. P - A) 

Tr[75(/^ + m) f'^ifi- A + H75(- P+ }i + m)] 

^ [fc2-m2 + i£][(A;- A)2-m2 + i£][(P-A;)2-m2 + i£]" ^ ' 

For circularly polarized (e+ = 0) initial and final photons* (/, J are t or J,), we obtain from Eq. (12) 

= (ei ■ )7+ + i{ei x el-^)37+75. (13) 

Equation (13) reduces to 7^(1 ±75) for the parallel helicities {i.e. + for tt and — for ||) and zero otherwise. Since 
the axial current 7"'" 75 does not contribute to the integral, i.e. ~ (kj_ x A_l) after the trace calculation, this term can 
be dropped. 

Thus, the DVCS amplitude {i.e. photon helicity amplitude) can be rewritten as the factorized form of hard and 
soft amplitudes 



M^^P,q,P') = M^^{P,q,P') = -el [ dx( + ^—)jr^{(,x,t), 

^ J \x — <^ + ie X — IS J 



(14) 



where 

iNc f dk~(fk± 



:FAC,x,t)= ='^J ^^^H^^Ak,P)K^Ak-A,P-A) 



Tr [75(/^ + m)7+(/^- A + m)75(- f + ^ + m)] 
[fc2-m2+ie][(A;-A)2-m2+ie][(P-A;)2-m2 + i£]" ' 



*As discussed in [16], for a longitudinally polarized initial photon, the Compton amplitude is of order 1/Q and thus vanishes 
in the limit — > 00. 
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The function J^^{(,x,t) is the "generahzed parton distribution" and it manifests characteristics of the ordi- 
nary(forward) quark distribution in the limits C ^ and t ^ 0. On the other hand, the first moment of J^jriC, x, t) is 
related to the form factor by the following sum rule [1,2]: 

/ dxJ'^{Q,x,t) = {l-Q/2)F^{t). (16) 
Jo 

In general, the polynomiality conditions for the moments of the GPDs [33,34] defined by 

dx x^-^TAi, X, t) = {l- C/2) i^„(C, t) (17) 

Jo 

require that the highest power of Q in the polynomial expression of Fn{C,t) should not be larger than n. These 
polynomiality conditions are fundamental properties of the GPDs which follow from the Lorentz invariance. The 
positivity property of the GPDs, following from positivity of the density matrix [35], has been discussed analyzing 
the double distributions (DD) [36]. While such conditions as sum rules and polynomiality arc more easily satisfied by 
the DD-based models for the GPDs, it has been noted [36] that the positivity constraints are more transparent in the 
framework based on the LF wavefunction. 

An important feature of the DVCS amplitude given by Eq. (14) is that it depends only on the skewncss parameter 
( = Q"^ /{2P ■ q) for large and fixed \t\{< Q^). This property of the exclusive DVCS amplitude is similar to the 
Bjorken scaling for the structure functions of the inclusive deep inelastic scattering (DIS) reaction. Moreover, for 
DVCS, the skewness parameter C, coincides with the Bjorken variable Xbj- Note also that, according to Eq. (14), the 
size of the imaginary part of the DVCS amplitude is proportional to T^r{(,^C,^t)^ to the value of the GPD at the 
crossover point x = C,. In its turn, the imaginary part of the DVCS amplitude determines the magnitude of the single 
spin asymmetry (SSA) [1] (see also [37,38]) that can be measured through the scattering of longitudinally polarized 
electrons on an unpolarized target. The SSA measurements for the proton target have been reported by HERMES 
[29] and CLAS [30] collaborations, with the magnitude of the SSA definetely inconsistent with zero. 

III. LIGHT FRONT BETHE-SALPETER APPROACH FOR THE GENERALIZED PARTON 

DISTRIBUTIONS. 

Let us now consider an explicit form of the GPD defined by Eq. (15) in the framework of Light-Front dynamics. 
Corresponding to the x < C, and x> C, regions of the covariant amplitude (see Fig. 3), the Cauchy integration over k~ 
in Eq. (15) has nonzero contribution coming either from the residue of the pole at {P— k)^ = rn? in the valence region 
a; > C or from the pole at k"^ = in the non-valence region x < C,. 




FIG. 3. Covariant amplitude with reduced photon vertex for pion GPD (left) and its non- valence x < C, part (right). 

These two cases correspond to two different time-ordered LF diagrams which we will consider separately. In the 

region of ^ < x the residue is taken at the pole in upper half-plane k~ relevant to (P — fc)^ = m?, i.e. k~ = 
P~ — (P — fc)~n + ie = P~ — [k^ + m?]/ (P — A;)+ -|- ie. Performing the Cauchy integration over k~ in this region we 
find 
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k^-m^ = k+k- - ki - ^ k+{p- - fc" - (P - A;)" ) = x (m^ - "^"^^ 

\ ( -L X j 

^2 _2 , A\+^D- A- /u n. A\-\ I ( k'l + vi^ 



{k - A)2 -m'^ik- A)+(P- - A- - (P - A:)-, - {k - A)-J = x' (^M^ - _ ^, J , (18) 
where 

a;' = and kV = k_L - A_l + x'A_i_ (19) 

are the "internal" LF-momcnta of the quark absorbed into the final pion state (x' evidently has the meaning of the 
plus-momentum of the absorbed quark measured in units of the final hadron plus- momentum). Thus, the Cauchy 
integration in Eq. (15) over k~ gives 



where 



^r'(C,^,i) = ^Y^(Y::^ / c!^k^x(2^2)(x,k^)5+,x'(2^2)(^',kl), 



2 , _2 u/2 , _2 



(20) 



X(2^2) - ^2 _ M^2 . ^^^0 - + ^ 

and is the trace in Eq. (15) evaluated at the point (P — fc)^ = m^; i.e. 

4P+ 

= 13^ (k^-k'^ + (21) 

The light-front vertex functions Hlf and h'j^p can be related to the covariant i?cov and -ff^ov vertex functions and 
then to LF wavefunctions, e.g., via instantaneous approximation following [18]. Here, following [20], we identify the 
LF wavefunction in X(2^2) as 



X(2.2)(x,kJ^^-^-L^^^^(.,k.), (22) 

where the Jacobian of the transformation k = {kz,'k±) (a;,k_L) is obtained using dkz/dx = Mo/[4a;(l — x)] with 
Mo = y x{i-x) ' ^^'^ radial wavefunction is given by 



<^(x,ki)=xj3^7^exp<j ^ ^ g^2 (23) 



The normalization factor x = 1 is introduced for future convenience and is fixed by the normalization of the pion 
form- factor F^{t = 0) = 1. After some simplifications, Eq. (22) can be written as 



X(2^2)(.,k,) = W — -^^(.,k,). (24) 



Substituting Eq. (22) into Eq. (20), we obtain the valence part of the pion GPD: 
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In the region of x < (, the residue shall be taken at the pole k"^ = rn? , i.e. k = 
located in the lower half of the complex-fc" plane. One can show that in this case 



[k^ + m^J/fc"*" — ze, which is 



{k - A)2 - 



x{l — x) 
^ '[ X"{1 - x") J 



where 



x" = x/C and k'j; = k_L + a;"A_L 



(26) 



(27) 



are the "internal" momenta of the quark annihilating into the photon. The second line in Eq. (26) is related to a 

"gauge-boson-wavefunction" [19]- 
Then the non-valence part of GPD after the Cauchy integration in Eq. (15) can be written as 



^™'(C,.,i) = ^1^ ,(1 / ^^k.X(...)(.,k.)5„+_,xn.,k",)/.',,, 



(28) 



where X2^2 and are the pion and gauge-boson LF wavefunctions and the final state LF vertex function h'^^p is now 
kept explicitly. The trace term at fc^ = w? {S^^^^^y) is given by 



4P+ 



(29) 




k 



P-k 



k-P 




FIG. 4. The difference in valence and non- valence LF vertex functions h'j^p and h'^p. 

Unlike in the valence part of GPD, the final state vertex h'j^p cannot be simply reduced to a LF wavefunction since 
it does not describe annihilation/creation of quark-antiquark pair out of the pion state. It corresponds to the decay of 
a quark into a quark and a pion state, or, equivalently the annihilation of a pion and a quark into a quark (see Fig. 4). 
However, it is worth noting that h'j^p for valence and non-valence contribution essentially are the special case of the 
same object, namely, the covariant vertex function H{k,P), in different kinematic regions. In particular, they obey 
the same BS equation and, in principle, can be deduced from the usual LF wavefunction by analytic continuation, as 
pointed out in Refs. [19,20]. 

Specifically, following Ref. [19], the x-functions in LF CQM are given by the solutions of the BS equation [24,39,40] 



{M^ - Ml)x{xi,ki±) 

= J [d'y][d^^±] /C(x,, k,^; yj, Ijj.) x(%, Iji), 



(30) 



where IC is the full BS kernel (which in principle includes all the higher Fock-statc contributions) and 
Mq = (to^ + k5_)/(l — x) — (m^ + k-]_)/x. Both valence and nonvalence BS amplitudes satisfy Eq. (30). For the usual 
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BS amplitude, referred as the valence wavefunction, x is greater than ^, while for the nonvalence BS amplitude x is 
less than C,. We use the notation for these two solutions 




FIG. 5. LF interpretation of the non-valence vertex (1 — > 3 blob). 

This notation is motivated by the relationship to the Fock state picture, in which the parton number before and 
after the kernel is interpreted for the nonvalence vertex as changing from 1 to 3. According to [19,20], the nonvalence 
BS amplitude shall be expressed via the valence BS amplitude and the full BS kernel in the relevant kinematic domain 
as illustrated in Fig. 5: 

= j M2/][c^^l-L]/C(a;i,ki_L;yj,lj_L)x(2^2)fe,lj_L)- (32) 

In a sense, the BS kernel in this region can be viewed as the sum of all processes resulting in creation of a quark- 
antiquark pair from the initial single quark state. With this in mind, we write the non-valence part of GPD (see left 
diagram in Fig. 6) as 



-^™'(C,.,*) = J^T^ / '^^k.X(...)(.,k,)5l„,,x^(.,k'V) 




FIG. 6. Non-valence part of GPD interpretation in LF CQM. 



It may be also noted that one should as well consider in this case the process in which quark-antiquark pair is 

created not by the spectator quark, but by the interacting quark (sec right diagram in Fig. 6). In Ref. [20], it was 
shown that this contribution can be reduced to the form of Eq. (33) with a simple modification of the kernel 



^(a;',k'i;j/',l^) = /C(ar',k'^;t/',l^) 



1 



'5nonval(?/'li)xi(?y,lj 



'^nonval('^' 

k_L)xL(a:,k±) 



(34) 
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where Xa,b denotes kinematic prefactors like x{l — x){l — x") etc. in the expressions corresponding to the diagrams 
shown in Fig. 6. 

Thus, we obtained the ampUtude corresponding to the nonvalence contribution given by Eq. (28) in terms of ordinary 
hght-front wavefunctions of hadron (x(2^2)) and gauge-boson (x^)- This method, however, requires the knowledge of 
the full BS kernel 1C{x, k_L ; y, which is in general dependent on the momenta connecting the one-body to three-body 
sectors, as depicted in Fig. 5. While the relevant operator K, in general involves all momenta (x, kj_; y, the integral 
of K, over y and 1_l in Eq. (33), 



dy' 



2/'(l-2/') 



d2l^^(a;,k^;y',l'^)X(2^2)(y',l'±), 



(35) 



depends only on x and k^. In this work, wc approximate G^r as a constant (mean value). This approximation has 
been previously tested in the analyses of exclusive semileptonic decay processes [19,20,24] and proved to be a good 
approximation at least in small momentum transfer region. The pion GPD, calculated in this way, looks like in Fig. 7. 



FIG. 7. The pion GPD at C = 0.5, t = -0.5 GeV^ computed using LFBS approach. 
Using the "effective" gauge-boson wavefunction in the non- valence region in the form 

one can ensure the continuity between the valence and the nonvalence parts of the GPD by suppressing the non-valence 
contribution at x' = x" = [20]. This cures the discontinuity found previously in Rcf. [19] whose origin can be traced 
down to the difference in the treatment of the valence and non-valence LF vertex functions. More specifically, the 
discontinuity in Ref. [19] originated from the crude approximation for the LF nonvalence vertex /i^^"™'. As discussed 
in the general form above, the nonvalence vertex /i^'™™' is an analytic continuation of h'^p into nonvalence kinematic 
domain so that it should be continuously and smoothly connected with the latter at the boundary between the valence 
and the nonvalence kinematic regions. We experimented with different forms of such extensions for the non-valence 
vertex, satisfying condition of continuity at x' = 0. One such possibility, e.g., could be 

(37) 

which has structure similar to that of the LF wavefunction, goes to a constant when x <^ and vanishes at x — > — 0, 
thus guaranteeing the continuity with the valence LF wavefunction hYp- However, the resultant behavior of the GPD 
in the non-valence region depends significantly on the detailed form of the model assumed for the extension of LF 
wavefunction into the non- valence region. On the other hand, while the introduction of the model for the nonvalence 



"8/3V'(l-a;") 



exp 



8^ 



(36) 



I, /nonval 



- m 



8a;/32 
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vertex satisfying these requirements is nontrivial, the inclusion of the virtual processes at the gauge-boson vertex using 
Eq. (36) removes the discontinuity at a; = C, thus alleviating the above problem at the lowest approximation of LFBS 
approach. Such a treatment ensures the vanishing of the nonvalence contribution at the crossover and also removes 
the infrared singularity in the amplitude. It was further suggested [20] that the value of the GPD at the crossover 
point needs not be zero because the higher Fock states in LFBS approach may introduce a nonzero contribution at 
this point. In the following section wo will consider this suggestion in more detail and show that, after inclusion of 
the qqg contribution into the DVCS amplitude, this value indeed is not zero. 



IV. EFFECT OF THE HIGHER FOCK STATES IN THE LFBS APPROACH FOR THE PION GPD. 

Our primary goal in this section is to analyze how GPDs are affected by the higher LF Fock states, especially at or 
near the crossover point x = C^. While the LFBS treatment in Section III can be considered as computing GPD with 
2-body Fock state contributions (upper diagram in Fig. 8), here we are interested in the effect of the 3-body Fock states 
as depicted in the lower diagram in Fig. 8. We concentrate specifically on the effect of the Fock states that include a 
gluon in addition to the constituent quarks in the DVCS amplitude. As suggested in Ref. [20], such contributions may 
result in a nonzero value at the crossover point of the GPD. 




FIG. 8. 2-body and 3-body contributions to pion GPD in LF. 

Given that there had been relatively little progress in finding exact, or even model, form of the 3-body LF wave- 
function in QCD from first principles, wc choose to model it by relating 3-body Fock states to 2-body Fock states. 
As shown in Fig. 9, we link 2-body states and 3-body states via a kernel including, in principle, all QCD processes 
which result in the production of one gluon along with the constituent quarks. Because of the complexity of such 
general kernel, in this work wc limit ourselves to model it by the lowest-order simplest possible process shown by the 
right diagrams in Fig. 9. Thus, when computing the GPD of the pion we will concentrate on the additional processes 
described by the covariant diagrams in Fig. 10. 
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FIG. 9. Relation between 2-body and 3-body Fock states (left) and an approximate relation for the kernel K (right). 




FIG. 10. Next-Fock-State corrections to DVCS which we denote SI, SI' and S2 (left to right). 



In Rcf. [18], it was stated that, when considering processes such as those of Fig. 10 in the LF dynamics, the time 
orderings when the additional gluon was exchanged entirely before (after) the photon emission should be omitted. 
This is because such contributions are related to initial (final) state interactions corresponding to iteration of the 
LF wavcfunction with BS kernel. Due to BS equation, these exchanges can be entirely absorbed back into the LF 
wavefunction. For example, the first diagram in Fig. 10, which we call SI, yields the three time-ordered contributions 
shown in Fig. 11. According to Ref. [18], two left diagrams in Fig. 11 are iterations of the LF wavefunction with the 
BS kernel in the one-gluon exchange approximation and, thus, should be absorbed into the LF vertex hYpiv, kj.). 




FIG. 11. Different time ordered contributions for a; > in SI. 



The only remaining contribution is the one where the gluon is present at the time of the photon emission, as shown 
by the right diagram in Fig.ll. Using the LF-perturbation theory, we write this contribution as 



/ 



IQ^3 167r3 -A; P-A) 



x(-C) X Tr[75(F P + m)-f"{^- p + m)j5{ji- A + m) + m)jail/+ m)] 

x(-l)x [l+k+{k-A) + {P-k)+{P-l)+{l-k)+]~'^ (38) 

{[P- - Ion -{P- l)on][P- - Kn " " Oon " " fc)o„] 
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[P- - A- - (fc - A)-„ - (P - Oon k)-J[P- - A- - (fc - A)- - (P - k)-J} ' . 

By the overline wc denote the parts which must be taken with the instantaneous contribution, i.e., with 
k~ = — X] This diagram contributes only to the region ( < x < y, and is zero otherwise. See Eq. (40) 
below for the definition of / and V. 

Similarly, for the final state interaction term SI', one can get three time-ordered diagrams in DGLAP region (modulo 
the instantaneous diagrams) shown in Fig. 12. According to Ref. [18], out of these three, only the left diagram in Fig. 12 
should be kept. Using the LF-perturbation theory rules, we obtain that this contribution is given by 



/ 



^^/^/..(.;P)/.o.(^-A;P-A) 



x(-CTr[75(^- P + m)ril/- P + m)j,{r- A + m)7a(^- A + m) ^{}i + m)]) 

X (-1) X [{I - A)+k+{k - A)+(P - fc)+(P - l)+{l - k)+] (39) 

[P- - A- - (fc - A)- - (P - Z)-„ - (/ - fc)-„][p- -A- -{I- A)-„ - (P - Z)- , 
where the overline has the same meaning as before. This term only contributes in the DGLAP region ( < x <y. 




FIG. 12. Different time ordered contributions for a; > f in SI'. 



In the ERBL region, there are two possible time-ordered processes (Fig. 13) which all should be retained as neither 
of them can be absorbed into the initial or final state vertex. The sum of these contributions yields the full covariant 
diagram, and we may use integration by poles, as described in greater details later, to compute this case for x < ( < y. 
Analogously, for the process described by the diagram 82 we keep all time-ordered contributions as, again, neither of 
them is absorbed into initial or final state wavefunction and we use integration by poles in this case as well. 




FIG. 13. Different time ordered contributions for x < ( in Si'. 

Combining these contributions, we calculate the GPD for ( = 0.3 and t = — 0.5GeV^. We use G^r = 0.32 and 
as 0.5. The results of our computation are presented in Fig. 14. As expected, !FTr{(,(^,t) ^ 0. However, we also 
notice a discontinuity in the GPD in Fig. 14 (the dashed line) near a; = C = 0.3. The origin of this discontinuity can 
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be attributed to the approximate nature of the model LF wavefunction that we used in our computation (namely, a 
Gaussian) while the Statement of Ref. [18] is only valid for the LF wavefunctions obtained from the exact solution of the 
BS equation. This shows that we cannot forget the time-ordered contributions omitted according to the prescription 
of Ref. [18] {i.e. the first two time-ordered diagrams in Fig. 11 and the last two diagrams in Fig. 12) but should 
include them all in our calculations. Taking into account these contributions can be seen as iterating our model LF 
wavefunction with the BS kernel once, to improve the quality of approximate wavefunctions. Of course, that would 
have no effect if we would already have the LF wavefunction as an exact solution of the BS equation in agreement 
with Ref. [18]. Since the covariant amplitude is expected to be continuous near x = and because retaining all 
time-ordered diagrams is equivalent to working with the covariant amplitudes, including these contributions is crucial 
to maintain the continuity of the GPD in the BSLF analysis with the model LF wavefunctions. 




FIG. 14. The pion GPD in LFBS approach with 2-body Fock states only (left), and with inclusion of specifically 
LF-time-ordered qqg contributions (right). The latter exhibits a discontinuity between DGLAP and ERBL region near x = 0.3 
(C = 0.3 and t = -0.5 GeV^). 



In order to consistently take into account the higher Fock states in our approach, we thus consider contributions 
from all processes in all time-ordered regions. We treat our calculations as corrections due to qqg Fock states in 
addition to the "effective" expression obtained in the LFBS approach of Ref. [19,20]. Then, each contribution in 
Fig. 10 corresponds to an expression in the covariant form 



where the overall coefficient is 



the effective vertex function is 



^ ^ / i2^r Denoimnitor [. . .] ' ^^^^ 



/ = Cpi-igsri^i-i) = -Cpg. 



V=li-^elN,){-g'^n 



C' " " ' \{q + k)^-m'^+is [k-q'f-m^+is 

and the trace is computed for the diagram with the effective- vertex contributing the 4 factor. We calculated the 
contribution from these processes using the LF momentum variables. First, we carried out integration by poles in 
Eq. (40). In the Cauchy integration over k~ and l~ , we used the poles located in the opposite halves of the complex 
k~ and l~ planes so that always a factor (27ri) or {—2-Ki) was introduced. After integration by poles, we are left with 
the expression 

C/ J 16^4 I67r4^ —V hLFh'LF(-Tr [...])/ (Denominator [.. .]) (41) 

1 

= N.i-ieDi-g^f J + ^) / '^-'^^ (-^^) / (Denominator!...]) , 
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where the overhne means that the expressions are taken with k~ and /~ values corresponding to specified poles for 
k~ and l~ in the complex plane. Comparing with iM'^'^ , iM^^ in Eq. (14) and noting that P~^Q^/Cq~ ~ 1 for DVCS, 
we identify the corresponding contribution to J^t^ from each process as 

fNe J ^ hLFh'i^p / (Denominator!...]) . (42) 

In this form, the tree-level contribution is simply 

hLph'^F (-C^H^/Q') / (Denominator [. . .]) , (43) 

and the higher Fock states corrections carry additional factor of f/lGn^ relative to the tree- level. The ± sign in 
Eq. (43) refers to two possibilities for the pole selection in the upper or lower halves of the complex plane, as we 
discussed in Section II. Each LF contribution can be constructed using Eq. (42) and the pole assignment for k~ and 
l~ presented in Appendix A. 

For each particular diagram in Fig. 10, wc obtain the following expressions. 
(SO) The diagram SO in Fig. 3 corresponds to the covariant expression 

* N.V ,.J-^?. H'^- ''.T^-'f- f„r"iy -^""l . , X HUt, P)H,..(t - A, P - A) , (44) 



(27r)4 (fc2-m2 + ie)((P-fc)2-TO2 + i£)((/!c-A)2-m2 + ie) 

so that the LF-exprcssion is given by 

, f dk {-C)TT[j5i^~ p + mh,ij^- |^ + m) fii^ + m)] 

y 167r3 " {k^ - + ie){{P - k f - + ie){{k - Ay - + is) ^''^'" ' ^ ' 

where k~ is set to its value at the corresponding pole. This case was in detail considered in the previous section. The 
± sign in Eq. (45) shall be read as plus for a; > C and minus for x < (. 

(SI) The initial-state-interaction contribution SI in the left diagram of Fig. 10 corresponds to the following expression 

r!k rll 

fVHi^{l;P)Hout{k-A;P~A) (46) 



(27r)4 (27r)4 

(-0 Tr [75(F yP + m)7"(/^- f + m)j5{fi- ^ + m) + m)7c.(^+ m)] 

^ ((P - 0^ - m2 + ie)((P - A:)2 - + ie){{k - A)^ - + ie){k'^ -m? + ie){P - tv? + ie){{l - kf + ie) 

Here, as well, k~ and l~ should be taken at the corresponding pole as presented in Appendix A. 
(SI') The final-state-interaction diagram SI' in the center diagram of Fig. 10 corresponds to 



/ 



dk dl 

(2^0^ (2^0^^^^'"^'^' -A-P-A) (47) 

(-0 Tr [75(^- f + m)7°(/- f + m)75(F 4 + m)7o(^- ^ + m) ^(^ + m)] 

^ ((P - kf - m2 + ie){{P - If - m2 + ie){{l - Af - + ie){{k - Af - m2 + is){k^ - + ie){{l - kf + ie) ' 

with the specific values for k^ and obtained from the poles presented in Appendix A. 
(S2) The box-diagram S2 in the right diagram of Fig. 10 is given by 

dk dl 

fVHinil; P)Hout{l - A; P - A) (48) 



(27r)4 (27r)4 

i-QTr [75(r- /> + m)75(r- A + m)7"(^- A + m) gf(^ + m)7^(/^+ m)] 

((P - Z)2 - m2 + ie){P -m?+ is){{l - Af -m? + is){{k - Af -m?+ ie){k^ - m? + ie){{l - kf + is) ' 

with the specific values for k~ and l~ obtained from the poles listed in Appendix A. Given the complexity of the LF 
expressions obtained after relevant poles substitution, we do not find it possible to present them more explicitly. The 
calculations for the pion GPD, including large portion of symbolic math, were further carried out numerically. 
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Now, let us discuss our numerical results for the pion GPD with the qqg contributions included. We performed our 
calculations for C = 0.1,0.3,0.5 and C = 0.7 and the value of -t = 0.5 GeV^. We used w 0.5 and = 0.32 [20]. 
While Gtt may be thought to depend on the values of ( in order to accommodate the sum rule for the GPD 

J dxT{C,x,t) = {l-CmF^{t), (49) 

we find, similar to Ref. [20], that this condition is satisfied within 10% for all considered values of C, even with G^r 
kept as a constant G^r = 0.32. Also, due to the additional contributions from 2-body and 3-body Fock states and the 
necessity to satisfy the form-factor normalization condition 

j dxJ^{C = 0,x,t = 0) = l, (50) 

we need to adjust the normalization factor for our model 2-body wavefunctions introduced in Eq. (23). We find that 
a factor of x « 1.5 is needed for the 2-body normalization to comply with Eq. (50) after qqg corrections are taken into 
account. 




FIG. 15. The pion GPD in LFBS approach including qqg corrections for —t = 0.5 GeV^ and different values of In the 
upper left, ( = 0.1, while C = 0.3 in the upper right, ^ = 0.5 in the lower left, and C = 0.7 in the lower right panel, respectively. 
The thick dashed lino represents the GPD calculation in the LFBS approach with only 2-body Fock states contribution, and 
the continuous line is the result for GPD after the qqg corrections were taken into account. 

As expected, due to the higlic;! Fock state contributions the value of the GPD at the crossover x = C is nonzero (see 
Fig. 15). Also, we now find that the GPD is continuous over the entire range of x and ( including x w C. It is crucial 
to take into account all possible time-ordered contributions including those that could be formally absorbed into the 
BS amplitude of the initial (final) vertex. After this was done, the connection between DGLAP and ERBL region is 
now continuous. 

We observed that for the smaller ( the higher Fock states introduce dominantly negative corrections in ERBL region 
so that J^iC, C, t) may become negative. Generally, we found that the GPD has sign-alternating structure, unlike the 
effective LFBS treatment result [20] and simple DD-based models [35,41], with nonvalence contribution x < ( being 
often negative and valence contribution x > ( usually positive, similarly to the results of the model in Ref. [42]. 
The point J^{x,(,t) = in our calculations always gets shifted from x = toward the larger values of x. For the 
larger C, however, the structure of the GPD changes and the GPD becomes positive for all x with no point at which 
^{x, C, t) = 0. We also find that the nonvalence contribution gets suppressed for the larger values of C relative to the 
case of small (. The positive valence part of the GPD becomes the dominant contribution in this case. 



16 



V. CONCLUSION 



Wc have taken into account the higher Fock state (qqg) contributions to the pion GPD, and verified that the value 
of the GPD at the crossover point x = C is indeed nonzero J>(C, ^) 7^ 0- First, we followed the statement of Ref. [18] 
and found that, although J-'-n-iC, Ci t) ^ 0, a, discontinuity near x = C, occurs due to the approximate nature of our model 
LF wavefunction. The explanation is that the statement of Ref. [18] is valid only for the LF wavefunctions obtained 
from the exact solution of the BS equation. Thus, wc had to retain the first two time-ordered diagrams in Fig. 11 
and the last two diagrams in Fig. 12. Taking into account these contributions is equivalent to iterating our model 
wavefunction with the BS kernel, and it is crucial to maintain the continuity of the GPD in our LFQM analysis. We 
thus included all the possible time-ordered diagrams equivalent to the qqg Fock state contributions shown in Fig. 10. 
We carried out integration by poles using the pole assignment summarized in Appendix A and numerically computed 
the pion GPD including qqg contributions. From this calculation, we found that the GPD is continuous and the value 
at the crossover point x = ^ is nonzero as expected. The essential finding of the paper, namely, the link between 
the nonzero GPD at the crossover point and the higher Fock-state contributions is not specific to the pion case but 
applicable also for the proton as well as other bound states. 
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APPENDIX A: POLE LOCATIONS FOR THE QQG CONTRIBUTIONS 

(SI) For the initial-state-interaction diagram SI (left in Fig. 10), depending on relative magnitude of {x, y, Q there 
are six different kinematic domains. In each of them poles for Cauchy integration are chosen as follows. 



Region 


Pole locations 


y < X <C 


- m2 = 0, {P 


-kf 


- 


= 0) & {P -m' = 


0, [k - Af 


-m'^ = 


0) 


y <C <x 


{P - m2 = 0, (P 


-kf 


- vn? 


= 0) 








C <y <x 


{P -m^ = 0, [P 


-kf 


- vn? 


= 0) 








X <y <C 


{{P-lf-m'^ = 


0, e 


- m? 


= 0) 








X <C<y 


{{P-lf-m'' = 


0, e 


2 

— TO 


= 0) 








C <x <y 


{{P-lf-m'' = 


0, e 


2 

— TO 


= 0) & {{P-l f- 


w? = 0, (fc 


-Af - 


m2 = 0) 



(SI') The final-state-interaction diagram SI' (center in Fig. 10). In each kinematic domain the poles are chosen as 
follows. 



Region 



Pole locations 



y < a; < C 

y <C <x 

C <y <x 

X <y <C 

X <C <y 

C <x <y 



all poles are in one half-plane, integral is zero 
all poles arc in one half-plane, integral is zero 

n2 _ 



0) 



iil- 


- A)2 


2 

— TO 


= 0AP- 


-kf 


{{P 


-If 


2 

— TO 


= 0, fc2 _ 


2 

TO 


{{P 


-If 


2 

— TO 


= 0, P - 


2 

TO 


HP 


-If 


- m? 


= 0, fc2 _ 





TO 



0, e 



m 



0) & {{P-lf-m^=0, (k-Af 



0) 



TO 



0) 
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(S2) The box-diagram S2 (right in Fig. 10). Poles are chosen as fohows. 



Region 



Pole locations 



- m2 = 0, (fc - A)2 - m2 = 0) 
all poles are in one half-plane, integral is zero 
all poles are in one half-plane, integral is zero 
((P - 0' - = 0, P - m2 = 0) & {{I - A)2 
((P - If - m2 = 0, P - m2 = 0) 

((P - 0^ - = 0, fc^ - = 0) fc ((P - - = 0, (fc - A)^ - = 0) 



y < X < C 
y <C <x 
C,<y <x 
X <y <C 
X <C <y 
C <x <y 



■ m 



2 = 0,e-m^ = 0) 



APPENDIX B: ORGANIZATION OF NUMERICAL CALCULATIONS 

The calculation of corrections to DVCS originating from addition of the qqg Fock states has been performed with 
the help of Mathematica program. Each contribution was specified by its covariant expression for the trace obtained 
from the diagrams in Fig. 10, the list of the denominator factors entering into Eq. (40), and the assignment of poles for 
each of 6 possible kinematic domains. As a result, the expressions for full amplitudes were constructed. Essentially, 
for each of the diagrams the expression like the following was generated: 



(-C)Tr[...]/Q2 
Denommator[. . .J 



(Bl) 



where we used for Hlf 



Hlf = 7= — (p{x, k_L). 



(B2) 



The wavefunction (t>{x,'kj_) is given by Eq. (23). For the LF wavcfunction of the gauge-boson, wo used the expression 
that is slightly different from Eq. (36). Specifically, we dropped from the argument of the exponent in k"^ -|- m?. 
This resulted in a slight shift of the nonvalence contribution in the direction of the larger x compared to Ref. [19]. 
However, this effect was rather insignificant. 

The amplitudes obtained in this way need to be integrated in 5 dimensions, i.e. (P'k±dycPl±. This can be done 
using a Monte-Carlo(MC) algorithm. To improve the efficiency of the MC integration, we analyzed and subtracted 
possible singularities from the amplitudes. In particular, almost each amplitude carried an IR-singularity at y = x, 
1_L = kj_. While these singularities are integrable. 



/* / 



1 



\1±\<A 



A{y -x)- (U 



kx)2 



= finite, 



(B3) 



they can seriously degrade the efficiency of the multi-dimensional integration, 
contributions in the form 



To avoid this, we picked out such 



M = 



(-C)Tr[...]/Q2 
Denominator[. . . 



hLph'] 



LF 



1 



3/=x, l± 



A(y-x)-(U-k^)2 



(B4) 



and integrated them analytically over dy and d^l±. The remaining two integrations over d^'k± have been done 
numerically with the built-in algorithms in Mathematica. The result was stored as the "IR-contribution" . In Fig. 16, 
we present these amplitudes, corresponding to SO, SI and S2. From Fig. 16, one can see how the IR-singularity 
subtraction works in practice. 
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FIG. 16. Examples of amplitudes SO, SI, S2 for some x and C plotted as a function of y. Gray line is the original amplitude 
and dashed line is the IR-piece subtracted from it. IR piece is partially integrated analytically. The remaining amplitude (black) 
is integrated numerically in 5D using MC method. 

In the reduced box-diagram (right in Fig. 10), the struck parton momentum k produces an UV-divergence because 
this amplitudes falls off like 1/k^ as |k_L| — » oo. This results in a logarithmic UV-divergence. Note that, since the 
original box diagram is UV-finite, this UV-divergence is fictitious and is caused by the approximation for the effective 
vertex used herein, in which the struck parton momentum was neglected and its propagator was replaced with l/Q^, 
a procedure that is admissible as long as <^ Q^. We may handle this singularity by using a cutoff of the order of 
Q^. For k'^ > Q^, the struck parton propagator should behave as 1/fc^. In principle, one may consider renormalization 
of this effective vertex to remove the dependence on A. However, we did not do this here because this term had only 
a weak dependence on the cutoff In A^ w InQ^. In numerical calculations, we took this contribution in the form 



Denominator [. . . 
and did the integration in kj_ analytically: 



Ba 



|ki|<A 



(ki + m2) 



2^2 • 



(ki + m2)2 



TT <^ ln(l + AVm^) 



A2 



(m2 + A2) 



(B5) 



(B6) 



The remaining expression was integrated over dy(Pl± and stored as "UV-contribution" . Finally, after these parts were 
subtracted, the residual part of the amplitude was integrated numerically in 5 dimensions using MC method. 

For the convenience of using Monte-Carlo, the integration over y was further reduced. Normally, the integration 
over y would go over three regions, e.g. y<x<(, x<y<( and x<(<yiix<(. These three integrals were 
rescaled, e.g. y = x'r] in y<x<(,y = x+{( — x)r] in x < y < ( etc, so that < < 1 and the three contributions, 
corresponding to these different regions, could be added together. Their sum, as a function of r], was the cumulative 
amplitude to be actually integrated with 5D Monte-Carlo (see Fig.l7). 



-2Q0 
-400 
- 600 
-8Q0 

-looo 

-1200 



□ .2 



□ . 4 



□ . 6 



□ , 8 



FIG. 17. Example of cumulative amplitude as function of 77 which was integrated with 5D Monte-Carlo. 



The MC integration was organized in a series of bunches, each bunch containing Nj, ~ 100 integrations using 
Mi 10, 000 points. Each bunch is treated as a "measurement" of the integral. When the test-points are independent 



19 



in the MC-integration, each measurement (bunch) is independent from each other and the general principles of statistics 
can be applied. For each bunch, thus, we estimated the dispersion at and the mean of the measurement distribution. 
In this way, we estimated the value of the integral as well as the statistical error in each bunch from MC integrations. 
All bunches finally were combined with appropriate weight factors to produce the final estimate of the integral and 
MC-integration error. We tested our integration with different limits of integration in (k_L,lj_) to see if any accuracy 
is systematically lost due to finite bounds of integration in kj_, We observed that, except for rapid decrease of the 
MC-efficiency, no noticeable change occurred when the integration region was enlarged. Special attention was paid 
to a; = ^ point since many of the amplitudes become singular in this case. Although all of these singularities are 
removable, we explicitly found the analytical limit in each expression for a; — > C either from the left or from the right, 
and accommodated this in the final graphs as left (right) limit points at x = 

Thus, with our program we generated three different numbers for each point x, These were "IR-contribution" , 
"UV-contribution" and the residual integrated amplitude. These were added to produce the final result 

M{im.l=Mo + f{MlR + Mres+BAMuv), (B7) 

where A4q is the Oth-order amplitude, Mir is IR-piece, Aiuv is UV-piecc of the amplitude and Ba is the factor shown 
in Eq. (B6) with « Q"^, where UV-cutoff A 5GeV was typically used. Mres is the MC-integrated remaining 
contribution. These results for ( = 0.3 are plotted in Fig. 18. 




FIG. 18. The plot of all contributions to scale for ^ = 0.3. The final result, "IR-contribution", "UV-contribution" and the 
5D residual integral are denoted by boxes, diamonds, stars and triangles, respectively. 

One final note is in order. In the LFBS approach of Ref. [19], there is a free parameter G^r which can be adjusted 
to fit the pion form factor as 

[\xJ'^{C,x,t) = {l-C/2)F^{t). (B8) 
Jo 

In principle, may be allowed to depend on t and (. In that case one should solve for from Eq. (B8). To facilitate 

this process, wc note that G^r appears at most linearly from the non-valcncc final-state vertex in the final expressions. 
In our numerical procedure we explicitly separated G° and G^ parts of the amplitudes and computed separately M'j 
and so that 

Mf = Mj + G^M}. (B9) 

Then G^r can be easily found from Eq. (B8) by solving a linear equation. In our calculations, however, we did not 
solve exactly for G7r(i, C)) but used G,r = 0.32 following the work of Ref. [20] where the justification can be found. 
Similar to Ref. [20], we find that the condition (B8) is satisfied well for all values of C even with a constant G,r- 
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